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ABSTRACT 

We examine the rotation rates, sizes, and star formation (SF) efficiencies of a represen- 
tative population of simulated disc galaxies extracted from the Galaxies-Intergalactic 
Medium Interaction Calculation (gimic) suite of cosmological hydrodynamic simula- 
tions. These simulations include efficient, but energetically feasible supernova feed- 
back, but have not been tuned in any way to produce 'realistic' disc galaxies. Yet, 
they generate a large number of discs, without requiring extremely high resolution. 
Over the wide galaxy stellar mass range, 9.0 < log lo [M»(M )] < 10.5, the simulations 
reproduce the observed Tully-Fisher relation, the rotation curves of disc galaxies in 
bins of stellar mass, the mass-size relation of disc galaxies, the optical rotation to virial 
circular velocity ratio ('V^t/V™'), and the SF efficiencies of disc galaxies as inferred 
from stacked weak lensing and stacked satellite kinematics observations. They also 
reproduce the specific star formation rates of ~ L* galaxies but predict too low levels 
of star formation for low-mass galaxies, which is plausibly due to the finite resolution 
of the simulations. At higher stellar masses, log 10 [M*(M Q )] > 10.6, the simulated 
galaxies are too concentrated and have too high SF efficiencies. We conjecture that 
this shortcoming reflects the neglect of feedback from accreting supermassive black 
holes in these simulations. We conclude that it is possible to generate a representative 
population of disc galaxies that reproduces many of the observed trends of local disc 
galaxies using standard numerical hydrodynamic techniques and a plausible imple- 
mentation of the "subgrid" astrophysical processes thought to be relevant to galaxy 
formation. 

Key words: galaxies: evolution — galaxies: formation — galaxies: general — galaxies: 
haloes — galaxies: stellar content — galaxies: structure 



1 INTRODUCTION 

Understanding the formation of normal disc galaxies re- 
mains one of the key challenges in astrophysics today. 
A general physical picture has b e en in place for some 
time (e.g., iFall fc Efstathioul Il98d : IWhite fc Freiikl Il99ll : 
iMo et al.ll 19981 ). but attempts to form realistic disc galaxies, 
and to test the general framework in detail, using ab ini- 
tio cosmological hydrodynamic simulations, have had only 
limited success thus far. Specifically, the simulations have 
historically yielded galaxies with too large of a fraction of 

* E-mail: mccarthy@star.sr.bham.ac.uk 



thei r baryons locked u p in stars (the "overcooling problem" , 
e.g.. lKatz et~ai1ll996u . sizes th at are too small (the "angu - 
lar momentum problem", e.g., INavarro" fc Steinmet3l2000h . 
and rotation curves that decline with radius (rather than 
being approximately flat as observed) due to th e presence of 
an ex cessively large spheroidal component fe.g-. lAbadi et all 
2003). It is plausible that all three of these problems are re- 
lated. 

It has become increasingly clear from many different 
lines of evidence that feedback plays a very important role 
in regulating the star formation (SF) of galaxies and also 
in setting their internal structure. Modelling feedback from 
supernovae (SNe) in simulations of galaxy formation is chal- 
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lenging, as the injected thermal energy tends to be radiated 
away before it can have any hydr odynamical effect (e.g., 
iKatz et al.lll996l : iBalogh et al1l200ll ). This is likely due to 
the fact that simulations still do not have sufficient resolu- 
tion (nor all the physics necessary) to model the interstellar 
medium correctly. Cosmological simulations must therefore 
resort to subgrid recipes to try to overcome the overcool- 
ing p r oblem (see the discussions in iDalla Vecchia fe Schavel 
2008, 2012), which is plausibly the fundamental underly- 
ing cause of the other problems ment ioned above (e.g., 
ISales et al.ll2010l ; IScannapieco et ai]|20f ll ). Typically, the SN 
energy is injected kine tically rather than thermally (e.g., 
iNavarro fc White! 1 19931 ) and/or radiative cooling is sup- 
pressed by hand for a sh ort time in the i mmediate aftermath 
of energy injection (e.g.. | Gerritsenlll997l'). althoug h other so- 
lutions exist fe.g.. |Palla Vecchia fc Schavell2012r ). 

In the past few years there has been a good deal of 
progress on the formation of more realistic discs galax- 
ies in cosmological simulations. A number of groups have 
been carrying out very high resolution "zoomed" simula- 
tions of individual systems primaril y using a thermal SN 
feedback subgrid prescript i on (e . g.. lOkamoto et alJ 2005, 
2010l;lGovernato et al.ll2007l.l2010l.l2012l; IStinson et al.ll2oTfj|; 
Agertz et all l201ll: iGuedes et al.l l201ll; lAvila-Reese et all 



201 it iMaccio et all l2012l : iBrook et all 120121 : 



see also 

Piontek fc Steinmeta l201ll for zoomed simulations with a 
kinetic feedback model). Generally speaking, these stud- 
ies have found that higher numerical resolution, together 
with increased density thresholds for SF, result in much 
more efficient feedback which, in turn, allows the simula- 
tions to reproduce many of the properties of normal disc 
galaxies, including reasonable SF efficiencies, disc sizes, and 
approximately flat rotation curves. Along with these suc- 
cesses, however, there remain important shortcomings, such 
as the apparent inability of current simulations (and semi- 
analytic models) to match the observed specific star forma- 
tion rate J] of disc galaxies both locally and at high redshift 
(e.g., lAvila-Reese et al.ll201ll ; iDave et al.ll201ll ). 

While the ability of recent high resolution simulatons 
to reproduce a number of important observables is certainly 
encouraging, caution is warranted in applying the under- 
standing derived from them to the galaxy population in 
general. Firstly, the feedback implemented in many of these 
studies has been tuned to produce 'realistic' disc galaxies. 
However, given the uncertain (and possibly degenerate) na- 
ture of subgrid processes, it is dangerous to claim success 
on the basis of a handful of simulated galaxies. The de- 
gree of realism should be judged on the basis of examin- 
ing many different observables over a wide range of mass 
and redshift. Secondly, the systems that are chosen to be re- 
simulated at very high resolution have typically, although 
not always, been selected on the basis of their merger his- 
tory - the common lore being that quiescent merger his- 
tories give rise to disc galaxies. However, the link between 
merger history and galaxy morphology i s complicated (e.g.. 



Qkamoto et al. 20051; IZavala et all 120081 ; IScannapieco et all 



20091 ; Sales et al. 120121 ') and, in any case, the fact that the 



majority of isolated normal galaxies in the local Universe 



1 The specific star formation rate is defined as the ratio of the 
star formation rate to stellar mass of a galaxy. 



are spirals implies that no such selection based on merger 
history should be necessary. 

An alternative approach, therefore, is to try to sim- 
ulate large, representative populations of galaxies without 
tuning the feedback or explicitly selecting particular merger 
histories in order to produce 'realistic' disc galaxies at the 
present day. The downside of this approach is that the reso- 
lution must be reduced, and so too are the kinds of questions 
one can ask of the simulations. Also, there is obviously no 
guarantee that reasonable disc (or elliptical) galaxies will be 
produced in the simulation, as there has been no explicit 
tuning to get the 'right' answer. At present, it is compu- 
tationally feasible to simulate large populations (~ 10 ) of 
normal galaxies with ~ kpc scale resolution. This is the ap- 
proach we adopt in the present study. 

In particular, we use the Galaxies-Intergalactic Medium 
Interaction Calculation (gimic) suite of cosmological hy- 
drodynamical simulations, which were carried out by the 
Virgo Consortiu m and are described in detail bv lCrain et alJ 
l|2009h (see also lSchave et al.ll2010t ). The suite consists of a 
number of large, nearly spherical regions (~ 20/i _1 Mpc m 
radius) extracted from the Millennium Simulation and re- 
simulated at higher resolution with gas dynamics, including 
subgrid prescriptions for metal-dependent radiative cooling, 
SF, chemodynamics, and a kinetic implementation of SN 
feedback. 

The resolution of GIMIC is somewhat lower than in the 
studies noted above that have looked at the formation of 
individual disc galaxies. Our work is therefore complemen- 
tary to these very high-resolution studies in that, instead 
of aiming for maximum resolution on a single system (or a 
handful of systems), we elect to run lower resolution sim- 
ulations of a large population of galaxies, so that we may 
compare our results to large observational samples of galax- 
ies and place them in a wider context. We note that a num- 
ber of recent zoom simulation-based studies have implied 
that very high resolution is required in order to form reason- 
able disc galaxies at all, but the accuracy of this statement 
depends crucially on the nature of the implementation of 
important subgrid physics (particularly S N feedback) in the 
simulations. As we have shown elsewhere (ICrain et al. 2010; 
iFont et all l201ll ; iMcCarthv et all 120121 : ISales et alJ I2012I ). 
the GIMIC simulations have no difficulty in forming disc- 
dominated (or even bulgeless) galaxies even though they do 
not have extremely high resolution (in any case, we show in 
the Appendix that our results are robust to changes in the 
resolution). 

In the present study, we examine the relation between 
stellar mass and rotation velocity of disc galaxies (i.e., the 
'Tully-Fisher' relation), rotation curves in bins of stellar 
mass of the disc galaxies, the relation between stellar mass 
and size of disc galaxies, the SF efficiency of both disc and 
elliptical galaxies, the relation between stellar mass and op- 
tical rotation to virial circular velocity ratio (V op t/V v i T ), and 
the relation between stellar mass and specific star formation 
rate (sSFR) . We compare the simulations to a variety of re- 
cent local observations over a wide galaxy stellar mass range 
of 9.0 < log 10 [M»(M s )] < 11.5. Throughout we adjust the 
observational s tellar mas s meas urements so that they are ap- 
propriate for a lChabrieri (|2003h stellar initial mass function 
(IMF) (e.g., we add 0.05 dex [subtract 0.1 dex] in the case a 
Kroupa ['diet Salpeter'] IMF), which is what is adopted in 
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the GIMIC simulations. For stellar masses de rived using the 
M/L- colour scalings of iBell et al] (|200ot ) or lBell fc de Jongj 
l|200ll ). we also apply a small adjustmen t to the mas s using 
the expr ession given in the appendix of lLi fc White! (|200S ) 
(see also iDutton et alj l201ll ) so that the stellar masses are 
consistent with those derived from more accurate 'SE P fit- 
ting' of five-band SDSS data (jBlanton fc Roweij|2007h . 

The present study is organised as follows. In Section 2 
we present a description of the GIMIC simulations. In Section 
3 we present our main results on the Tully-Fisher relation, 
rotation curves, mass-size relation, and SF efficiencies. In 
Section 4 we present some scaling arguments that link the 
faint-end slope of the galaxy stellar mass function to the 
slope of the Tully-Fisher relation. We summarize and discuss 
our findings in Section 5. 



2 SIMULATIONS 

A detailed des c riptio n of the GIMIC simulations can be found 
in lCrain et al.l (|2009l ). We therefore present only a brief sum- 
mary here. 

The suite consists of a set of hydrodynamical re- 
simulations of five nearly spherical regions (~ 20/iT 1 Mpc 
in radius) e xtracted from the (d ark matter) Millennium 
Simulation l|Springel et alj 120051 ). The regions were se- 
lected to have overdensities at z = 1.5 that represent 
(+2, +1, 0, — 1, — 2)<x, where a is the root-mean-square de- 
viation from the mean on this spatial scale. The 5 spheres 
are therefore diverse in terms of the large-scale structure 
that is present within them. For the purposes of the present 
study we select 'central' galaxies only, which are defined as 
the most massive self-gravitating collection of stars, gas and 
dark matter (i.e., the most massive subhalo) in a friends- 
of-friends group. Systems and their substructures are iden- 
tified in the simu lations using the Subfind algorithm of 
iDolag et al] d2009l), th at ex tends the standard implementa- 
tion of lspringel et al.l (|200lT ) by including baryonic particles 
in the identification of self-bound substructures. 

The cosmological parameters adopted for GIMIC are the 
same as those assumed in the Millennium Simulation and 
correspond to a ACDM model with Q, m — 0.25, Qa = 0.75, 
Qb = 0.045, as = 0.9 (where erg is the rms amplitude of 
linearly evolved mass fluctuations on a scale of 8/i -1 Mpc 
at z = 0), H = WOh km s" 1 Mpc" 1 , h = 0.73, n s = 1 
(where n s is the spectral index of the primordial power 
spectrum). The value of erg is approximately 2-sigma higher 
than has been inferr ed from the most recent CMB data 
|Komatsu et al.l 1201 lh ■ which will affect the abundance of 
the simulated galaxies somewhat, but should not signifi- 
cantly affect their internal properties. 

The simulations were evolved from z = 127 to z — 
using the Tree PM-SPH code gadget-3 (last described in 
ISpringeli r2005). which has been modified to incorporate 
new baryonic physics. Radiative cooling rates for the gas 
are computed on an element-by-element ba sis using pre- 
computed tables generated with CLOUDY (jFerland et al.l 
1998) that contain cooling rates as a function of density, tem- 
perature, and redshift and that account for the presence of 
t he cosmic microwave back ground and photoionisation from 
a lHaardt fc Madau] (l200l l ) ionising UV/X-Ray background 
(see Wiersma et al.l2009a ). This background is switched on 



at z = 9 in the simulation, where it 'reionises' the simulation 
volume. 

Star formation is tracked in the simulat i ons f ollow- 
ing the prescription of ISchave fc Dalla Vecchial |2008). Gas 
with densities exceeding the critical density for the onset 
of the thermo-gravitational instability is expected to be 
multiphase a nd to form stars (uh ~ 10~ 2 — 10" 1 cm" 3 ; 
ISchave 1120041 ). Because the simulations lack the physics to 
model the cold interstellar gas phase, an effective equa- 
tion of state (EOS) is imposed with pressure P oc p 4//3 
for densit iefl n_g > n t where n« = 0.1 cm -3 . As de- 
scribed in ISchave fc Dalla Vecchiial (|2008l ). gas on the effec- 
tive EOS is allowed to form stars at a pressure-dependent 
rate that reprod uces the observed Kennicutt-Schmidt law 
|Kennicutdll998l) by construction. Note that the resolution 
of the GIMIC simulations (see below) was actually chosen 
so as to marginally resolve the Jeans mass at the adopted 
density threshold for SF. Higher resolution is necessary to 
be able to increase the SF threshold, but it is not suffi- 
cient by itself since additional physics is required to accu- 
rately capture the phase transition to the cold ISM (cer- 
tainly radiative transfer but possibly also dust, conduction, 
cosmic rays, and magnetic fields). The recent very high res- 
olution zoomed simulations discussed in Section 1, which 
have pushed to higher SF thresholds, have not included this 
additional physics. 

The timed release of individual elements by both mas- 
sive (Type II SNe and stellar winds) and intermediate mass 
stars (Type la SNe and asymptotic giant branch stars) is in - 
cluded following the prescription of IWiersma et all (|2009bl ). 
A set of 11 individual elements are followed in these sim- 
ulations (H, He, C, Ca, N, O, Ne, Mg, S, Si, Fe), which 
represent all the important species for computing radiative 
cooling rates, in the presence of a UV/X-ray background. 

Feedback from SNe is incorporated usin g the kinetic 
wind model of iDalla Vecchia fc Sc havc (2008) with the ini- 
tial wind velocity, v w , set to 600 km/s and the mass-loading 
parameter (i.e., the ratio of the mass of gas given a velocity 
kick to that turned into newly formed star particles), r), set 
to 4. This corresponds to using approximately 80% of the 
total energy available from SNe for a Chabrier IMF, which 
is assumed in the simulation. Note that unlike many previ- 
ous works that use similar feedback schemes, the "kicked" 
particles are not hydrodynamically decoupled at any stage. 
This leads to efficient entrainment of the ISM by outflows 
which can significantly affect the properties of galaxies (see 
Dall a Vecchia & Schavc 2008). The above choice of feedback 
parameters yields a good match to the peak of the SF rate 
history of the universe. The simulation feedback parameters 
were chosen without any regard to observables other than 
the peak of the cosmic SF history. 

To test for numerical convergence, the GIMIC simula- 
tions have been run at three levels of numerical resolution: 
'low', 'intermediate', and 'high'. The low-resolution simu- 
lations have the same mass resolution as the Millennium 



2 Gas particles are only placed on the EOS if their temperature 
is below 10 5 K when they cross the density threshold and if their 
density exceeds 57.7 times the cosmic mean. These criteria pre- 
vent SF in intracluste r gas and in intergalactic gas at very high 
redshift, respectively jSchave &: Dalla Vecc hia 2008). 
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Simulation, while the intermediate- and high-resolution sim- 
ulations have 8 and 64 times better mass resolution, respec- 
tively. Only the —2a and Oct volumes were run to z = at 
high resolution, owing to the computational expense of these 
calculations. We adopt the high-resolution simulations in our 
main analysis and make comparisons to the intermediate- 
resolution simulations to test for numerical convergence. 
The high-resolution runs have a dark matter particle mass 
7?idm — 6.62 x 10 6 /i _1 Mq and an initial gas particle mas^f] 
of m gas ~ 1.45 x 10 6 /i -1 Mq. The Plummer equivalent grav- 
itational softening is 0.5/i _1 kpc, which is fixed in physical 
space for z ^ 3. 

We present a convergence study in the Appendix which 
shows that the simulation results pertaining to star for- 
mation efficiency, kinematics, and sizes are converged for 
systems with at least ~ 100 star particles. In the analy- 
sis below, we adopt a conservative lower stellar mass cut of 
log lo [M*(M0)] = 9.0, which corresponds to approximately 
500 star particles. We note, however, that a larger number of 
star particles appears to be required before converged sSFRs 
are obtained (see the Appendix). 

Lastly, the GIMIC simulations are based on the SPH for- 
malism. A relevant issue is whether differences in simulation 
techniques (e.g., SPH, fixed mesh, moving mesh) are impor- 
tant. It has been shown using idealised tests that standard 
SPH (which was adopted for the GIMIC simulations) has a 
number of problems relating to resolvin g fluid instabilities 



and the treatment of weak shocks (e.g., Agertz et alj 20071; 
Mitchell et al] 120091 ; IVogelsberger et all 12011 " ISiiacki et al.l 
20121 ). While it would be prudent to address these issues in 



the future (perhaps using modifications to the SPH fo rmal- 
ism that appear to mostly r esolve these issues; e.g., iPricd 
120121 ; iRead fc Havfieldll2012|), we point out t hat both the 
OWLS llSchave et all I2OI0I : ISales et ail l201Ch and Aquila 
|Scannapieco et al* l201ll )~proiects have clearly shown that 
the differences obtained using different hydrodynamic tech- 
niques are small compared to the differences in the results 
that arise from the use of different subgrid prescriptions. 



3 RESULTS 

Below we examine the z = stellar mass-rotation veloc- 
ity ('Tully-Fisher') relation (Section 3.1), rotation curves 
in bins of stellar mass (Section 3.2), the mass-size relation 
(Section 3.3), the SF efficiency of the simulated galaxies 
(Section 3.4), the stellar mass-V op t/K,i r relation (Section 
3.5), and the stellar mass-sSFR relation (Section 3.6). In 
each case we compare our results to a variety of recent lo- 
cal observations. We focus on the wide stellar mass range, 
9.0 ^ log lo [M*(M0)] ^ 11.5 (which corresponds roughly to 
a halo mass range of 11.3 < log 10 [M2oo(Mo)] < 12.7 in the 
simulations, where M200 is the total mass within the radius, 
i~200, that encloses a mean density of 200 times the criti- 
cal density). The lower stellar mass limit is motivated by 
the numerical convergence study in the Appendix, while the 
upper mass limit corresponds to the most massive galaxies 



3 When a star particle forms, it has the same mass as the gas 
particle from which it formed (the gas particle is removed). Thus 
the initial stellar mass is typically ~ 10 6 h~ 1 Mq. However, star 
particles lose mass over cosmic time due to stellar evolution. 



in the simulation volumes. Within this mass range, there are 
714 simulated galaxies in total in the GIMIC ( — 2,0)<r high- 
resolution volumes. This includes both disc- and spheroid- 
dominated galaxies. 

When making comparisons to the observed Tully-Fisher 
relation, the mass-size relation, the observed rotation curves 
of disc galaxies, the V op t/V v i T , and the sSFRs, it is appropri- 
ate to select only disc-dominated systems from the simula- 
tions. One way to do this is kinematically; e.g., to use a cut 
in the fractional contribution of ordered rotatio n to the total 
kinetic energy and/or angular momentum (e.g-. lAbadi et all 
l2003l ; [Sales et all2010l ). However, in general it is not possible 
to do this kind of selection for observed galaxies. Therefore, 
we instead adopt a simple procedure for assigning galaxies 
to either the disc- or spheroid-dominated categories based 
on the distribution of their stellar light. For each galaxy we 
produce a projected azimuthally-avera ged surface brigh tness 
profil<0in the SDSS i-band (e.g., as in lshen et ai1l2003l ). We 
fit this radial profile with a Sersic law: 



fl(R) = fl e + 



2.5b n 
ln(10) 



{R/R t 



(1) 



where R is the projected distance from the galaxy center, 
R e is the effective radius, /j e is the surface brightness at 
R e , n is the Sersic index, and b n is a func tion of n only 
that must be computed numerically (see, e.g-. lGraham et all 
2005). Finally, to approximately mimic the depth of SDSS 
observations, we restrict the fit to radii where jJ,(R) < 24.0 
mags, arcsec -2 . The results are not sensitive to the exact 
surface brightness cut adopted. 

We use the best-fit Sersic index, n, to assign galaxies 
to the disc-dominated and spheroid-dominated categories. 
The case of n = 1 corresponds to a purely exponential dis- 
tribution, which is known to describe the disc component 
of galaxies very well, while n = 4 corresponds to the classi- 
cal de Vaucouleurs distribution, which has historically been 
used to describe the surface brightness profiles of e lliptical 
galaxies. It has been shown bv lBlanton et all l|2003bl ) using a 
large sample of galaxies in the SDSS that there is an approx- 
imately bimodal distribution in n at fixed galaxy luminosity 
(or stellar mass), in analogy to the well-known bimodal dis- 
tribution of galaxy colours at fixed luminosity. Red galaxies 
are typically those with large Sersic indices (n > 3), while 
blue galaxies correspond to those wi th relatively small val- 
ues of n < 2 on average. Following IShen et all (|2003l ). we 
use n = 2.5 as our cut for distinguishing between disc- and 
spheroid-dominated systems. With this cut, there are 432 
disc-dominated and 282 spheroid-dominated galaxies over 
the range 9.0 < log lo [M»(M )] < 11.5 (see Fig. H below, 
which shows the distribution of n with galaxy stellar mass) . 



4 The i-band luminosities are computed by treating each star par- 
ticle as a simple stellar population (SSP). The simulations adopt 
a Chabrier IMF and store the age and metallicity of the particles. 
We use these to compute a spectral energy d istribution for each 
star p article using the GALAXEV model of iBruzual fc Charlotl 
( 2003). The i-band luminosity is obtained by integrating the prod- 
uct of the SED with the SDSS i-band transmission filter function. 
Our calculation neglects the effects of dust attenuation. 
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3.1 Stellar mass - rotation velocity relation 

In Fig. [TJ we examine the relatio n betwe en galaxy stellar 
mass and circular velocity (V c = y GM/r) for disc galaxies, 
which is commonly referred to as the 'Tully-Fisher' (here- 
after TF) relation. The top panel shows the relation for the 
case where the circular velocity is measured at the 3D radius 
that contains 80% of the total i-band flu43 in projection. We 
refer to this velocity as Vso- The bottom panel shows the re- 
lation between the galaxy stellar mass and the maximum 
circular velocity, Knax- 

We focus first on the top panel of Fig. [TJ The cyan points 
with error bars co rrespond to the recent measurements of 
iReves et all |201lf ). These authors selected a broad sample 
of galaxies from the SDSS and derived rotation curves from 
follow-up Ha observations (so me of whic h were already pre - 
sented in lPizagno et al.ll2007l ). Although IReves et alj (20111 ) 
did not focus specifically on disc galaxies, the requirement 
that the galaxies have usable Ha rotation curves effectively 
weeds ellipticals out of their sample. We use their estimates 
of the rotation velocity, Vgo, and their published r-band 
magnitudes and g — r colours t o convert t o stel lar mass 
using the M/L colour scalings of Bell et al.l (|2003l ). We se- 
lect only galaxies from IReves et al.l ( 2011 ) with axis ratios 
0.3 < b/a < 0.5 in order to minimise inclination uncerta in- 
ties and extinction corrections (as in lDutton et al.ll2010f) . 

Over the wide range 9.0 < log 10 [Af» (M©)] < 10.5 the 
simulated galaxies fall approximately on top of the observed 
relation. To our knowledge, this is the first time a cosmo- 
logical hydrodynamic simulation has produced a population 
of galaxies that is consistent with the observed TF rela- 
tion. The agreement over this mass range is not perfect 
(the observed relation appears to be slightly steeper than 
the simulated one), but it is clearly a very large improve- 
ment over previous simulat ion studies that looked at r epre- 
sentative populations (e- g., Navarro fc Steinmet j|200ol ). We 
note that Ide Rossi etail (|20ld) have also recently examined 
the TF relation of a population of galaxies simulated with 
slightly higher resolution but in smaller boxes (10/i -1 Mpc 
on a side). Their simulated TF relation has approximately 
the correct slope but no direct comparison was made with 
observations so it is unclear whether the simulation repro- 
duces the observed zero point as well. 

At higher masses (log 10 [M*(Mo)] > 10.6), the simu- 
lated galaxies have circular velocities that are significantly 
larger than the observed rotation velocities. This is likely 
the result of overcooling in massive halos (indeed, see Fig. [4] 
below). This may signal that an increase in the efficiency of 
SN feedback is required at high masses and/or that another 
form of feedback (e.g., from Active Galactic Nuclei, here- 
after AGN) is energetica lly important at th ese mass scales 
(e^^ 



Benson et al.ll2003l ). On the other hand. Ide Rossi et al.l 



(2010), who also implemented only SN feedback, did not find 
such a sharp change in the slope of their TF relation at high 
masses, although due to their smaller simulation box size 
they had a much smaller number of high mass galaxies in 
their sample. In practice, the mass scale at which overcool- 
ing begins to be an issue will depend on the details of the 
subgrid prescriptions and the adopted parameters, as well 

5 We use the best-fit Sersic model to determine the projected 
radius that encloses 80% of the flux. 
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Figure 1. Comparison of simulated and observed relations be- 
tween galaxy stellar mass and rotation velocity for disc galaxies 
(i.e., the 'Tully-Fisher' relation) at z = 0. The top panel shows 
the case where the velocity is measured at the radius that encloses 
80% of the total i-band flux. The bottom panel uses the maximum 
velocity. Circular velocities are used for the simulated galaxies. 
The simulated galaxies lie approximately on top of the observed 
relation over the wide range 9.0 < log 1( )[Af, (Mq)] < 10.5. At 
higher masses, the simulated galaxies have circular velocities that 
are too large, which we show later to be due to overcooling in the 
largest systems. 



as on the resolution. However, as many studies have demon- 
strated, it is not obvious that the problems at the massive 
end can be solved with only feedback from star formation. 

It also appears that the simulated TF relation has 
smaller scatter than the observed one. However, it is not 
clear whether the difference in scatter is a genuine discrep- 
ancy: since we have measured circular velocities, rather than 
gas rotation velocities (which will not be perfectly circular) , 
the scatter in the simulated TF relation is likely to be un- 
derestimated. Particle noise (due to our finite resolution) 
inhibits the accuracy to which the gas rotation velocities 
can be determined in the simulations, which is why we have 
used circular velocities instead. In addition, scatter in the 
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observed TF relation can be introduced by errors in correc- 
tions for disc inclination and dust attenuation (although we 
have attempted to mitigate this by selecting observed galax- 
ies within a particular axis ratio range), as well as errors in 
distance measurements. A much more careful analysis, tak- 
ing into account galaxy selection criteria and observational 
uncertainties, is required to make quantitative statements 
about the degree of similarity (or lack thereof) in the scat- 
ter in the observed and simulated galaxy populations. 

In the bottom panel of Fig. [1] we compare the observed 
and simulated TF relations using Vmax rather than Vso- 
The red poi nts with error ba rs represent the 21 cm mea- 
surements of IVerhel icn (2001), who measured the TF rela- 
tion for a volume-limited complete sample of spiral galax- 
ies in the nearby Ursa Major Cluster. We derive stellar 
masses for their galaxies using their /-band magnitudes and 
B—R colours and the M/L colour scalings of lBell fc de Jong 
d200ll). We als o show the best-fit power-law relatior0 of 
iBell fc de Jond (|200lh base d on K-band data (dash ed or- 
ange line), as well as that of lAvila-Reese et alj l|2008l ) based 
on a compilation from the literature (dot-dashed green line 
represents their 'orthogonal' best fit and dotted green lines 
represent ±la intrinsic scatter). 

In accordance with the results shown in the top panel 
of Fig. [T] we find that the simulated TF relation (using 
Vinax) agrees well with the observed one over the wide range 
9.0 < log 10 [Af„(Mo)] < 10.5 while at higher masses the sim- 
ulated galaxies have too high circular velocities with respect 
to observed disc galaxies of the same mass. 

In terms of the degree of scatter in the V ma x — Af* re- 
lation, we note that the role of distance un certainties in the 
observed relation should be minimal for the Vcrhciicn (2001) 
results, since their galaxies belong solely to the Ursa Major 
cl uster. Inter e stingl y, the scatter in the Vma x — M * relation 
of IVerheiienl |200ll ) and lAvila-Reese et al.1 |2008l ) appears 
to be similar to that of the simulations for galaxies with 
Iog lo [M*(M )] < 10.0. 

A caveat to bear in mind when comparing the simu- 
lated and observed Vmax — M, relations, is that Vmax for the 
simulated galaxies is the maximum circular velocity within 
»"200, whereas for the observations it is the maximum rota- 
tion velocity within the region out to which it is possible 
to measure gas rotation speeds (typically ~ 20 kpc). We 
have investigated what happens to the simulated V ma x-M* 
relation when we limit the maximum radius to 2 effective 
radii (to crudely mimic the extent of the HI gas). We find 
a nearly identical relationship for log lo [Af»(M0)] > 9.5 (the 
measured Vmax is the true one). However, at lower masses, 
the radius where the circular velocity curve peaks is beyond 
2 effective radii and, consequently, the measured Vmax is 
lower than the true one, but only by « —0.05 dex on aver- 
age. This has the effect of shifting the low mass simulated 
galaxies down so that they tend to li e in between th e re- 
sults of lAvila-Reese et all (|2008l ) (and IVerheiienll200ll ) and 

6 Strictly speaking, this refers to the Vfl a t — Vf* relation of 
I Bell fc de Jonel 112001! ), where Vfl a t is the best-fit asymptotic ve- 
locity at large radii. For galaxies with flat rotation curves Vflat 
is obviously equivalent to V ma x. For galaxies with rising rotation 
curves, e.g., low-mass galaxies, Vfl a t > V m ax as Vmax is deter- 
mined from the observed rotation curve whereas Vfl a t represents 
an extrapolated quantity. 



IBell fc de Jond (|200ll ). yielding even better agreement be- 
tween the simulations and the observations. 

As discussed at the beginning of Section 3, we selected 
the disc-dominated systems from the simulations for com- 
parison to the observed TF relation (of disc galaxies). How- 
ever, it is worthwhile to comment briefly on the behaviour 
of the spheroid-dominated galaxies, which are not shown 
in Fig. [T] Over the mass range where the disc-dominated 
simulated galaxies agree well with the observations (9.0 < 
log 10 [M*(Mo)] < 10.5), the spheroid-dominated galaxies 
follow the disc-dominated trend closely. In other words, our 
results are insensitive to the cut in Sersic index for this mass 
range. It is only at the highest masses where we see a no- 
ticeable difference, in the sense that the spheroid-dominated 
galaxies have lower circular velocities (Vso and V m ax) at a 
fixed stellar mass. However, as we will show in Section 3.3, 
it is clear that both disc- and spheroid-dominated galaxies 
at high stellar masses suffer from strong overcooling with 
respect to observed galaxies. 

The TF relation encapsulates the rotation of a disc 
galaxy with a single number. However, since neither sim- 
ulated nor observed galaxies have perfectly flat rotation 
curves, it is possible for simulations to lie on the observed TF 
relation without reproducing the rotation curves of real disc 
galaxies. Furthermore, it has been argued that the shapes of 
the rotation curves of normal disc galaxies ar e inconsistent 
with the predictions of the ACDM model (e.g. JSalucci et all 
2012), as deduced from dark matter only simulations. With 
this motivation, we next examine the rotation curves of ob- 
served and simulated galaxies as a function of galaxy stellar 
mass. 



3.2 Rotation curves in bins of stellar mass 

In Fig. [2] we compare observed galaxy rotation curves with 
predicted galaxy circular velocity curves in 9 equal-width 
bins of log 10 M* spanning the range 9.0 ^ log 10 [M*(Me)] ^ 
11.25. The continuous black curves (connecting filled black 
circles) represent the median circular velocity profiles of the 
simulated galaxies and the error bars represent the 5 th and 
95 th percentiles. The dot-dashed red, dashed purple, and 
dotted orange curves represent, respectively, the median cir- 
cular velocity profiles of the dark matter, stars, and gas 
separately. The shaded vertical region delineates the range 
^ r ^ 3 softening lengths (~ 2.05 kpc), where the circular 
velocities of the simulated galaxies are expected to be un- 
reliable. The cyan cu r ves re present fits to the Ha rotation 
curves of I Reyes et al.l (|201ll ). which we reproduce using th e 
best-fit arctan parameters provided by iReves et al.l (|201lf ). 
The best-fit curves are plotted out to the radius that en- 
closes 80% of the i-band flux. 

The agreement with the amplitude and shape of 
the observed rotation curves over the mass range 9.0 < 
log 10 [M* (M Q )] < 10.5 is encouraging. The observed and 
simulated velocity curves both show a gentle rise from ~ 2 
kpc out to about 5 kpc, on average, before levelling off to an 
approximately constant value at larger radii. We note that 
over this range in galaxy stellar mass, dark matter is the 
most important contributor to the circular velocity curves 
in the simulated galaxies (except in the innermost regions of 
the last mass bin), as can easily be seen by comparing the 
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Figure 2. Observed rotation curves and predicted galaxy circular velocity curves in bins of galaxy stellar mass. The continuous black 
curves (with filled black circles) represent the median circular velocity of the simulated galaxies and the error bars indicate the 5 th and 
95 th percentiles. The dot-dashed red, dashed purple, and dotted orange curves represent, respectively, the median circular velocity profiles 
of the dark matter, stars, and gas separately. The shaded vertical region delineates the range ^ r ^ 3 softening lengths (f» 2.05 kpc), 
where th e circular ve locities of the simulated galaxies are expected to be unreliable. The cyan curves represent fits to the Ha rotation 
curves of lReves et al 1 l|201lh . Over the range 9.0 < log 10 [M»(MQ)] < 10.5 the predicted circular velocity curves are established primarily 
by the dark matter and look remarkably similar to the observed rotation curves. The simulated galaxies with log 10 [M* (Mq )] > 10.6 are 
rotating too rapidly in the inner regions, where the stellar component dominates the dynamics, due to overcooling (see Fig. [!}. 



dot-dashed red curve (dark matter) to the dashed purple 
(stars) and dotted orange (gas) curves. 

Strong deviations of the simulated circular velocity 
curves from the observed rotation curves are evident at 
higher stellar masses, however, in the sense that the sim- 
ulated galaxies have mass distributions which are too con- 
centrated, particularly in the very innermost regions (r < 5 
kpc). This is clearly due to the stellar distribution being 
much too compact and dominant in these systems. It is 
worth pointing out that many observed massive galaxies do 
show evidence fo r declining rota t ion velocity cu rves (e.g., 
iPersic et al. 1996]; ICourteau|[l997l ; IVerheiienll200lh and that 
such behaviour cannot be captured by the arctan functional 



form adopted bv lReves et alJ (|201ll ) in fitting their rotation 
curves. However, it is very unlikely that this accounts for 
most of the difference at high stellar masses. Indeed, the 
fact that for these high-mass galaxies there is a normalisa- 
tion offset at relatively large radii, strongly suggests that 
these galaxies have too high SF efficiencies. 

A comparison of the rotation curve results to the TF 
results in Section 3.1 reveals a nice consistency, in that the 
rotation curves look sensible over the range of stellar mass 
for which the simulated galaxies lie along the TF relation, 
and vice versa. 

In terms of the spheroid-dominated galaxies in our over- 
all sample (which are not used in Fig. [2} , there are no sig- 
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nificant differences between disc-dominated and spheroid- 
dominated galaxies in the circular velocity curves until high 
stellar masses (M* > lO n M0) are reached. At high stel- 
lar masses the circular velocity curves of the spheroid- 
dominated galaxies are systematically lower in amplitude 
and somewhat flatter than for the disc-dominated galaxies 
of the same stellar masfl 

3.3 Stellar mass - size relation of disc galaxies 

In Fig. we plot the relation between galaxy stellar mass 
and size. We use the radii which enclose 80% (top panel) and 
50% (bottom panel) of the i-band flux in projection. The 
filled black circles with error bars correspond to the median 
and 16 th and 84 th percentiles in bins of galaxy stellar mass 
for the simulated GIMIC galaxies. The radii R50 and Rgo are 
derived by fitting a Sersic law to the azimuthally-averaged 
i-band surface brightness profiles of the simulated galaxies, 
as described above. For comparison, we plot the median 
M« — -R50 disc relations (and 1-sigma s catter) of IShen et al.l 
|2003l ) (derived from SDSS data) and IBaldrv et al l (j2012T ) 
(derived from GAMA data) in the bottom pa nel, and in- 
dividu al measurements of M* and Rgo from iReves et al.l 
ll201ll) in t h e top pane l. Note that t h e com parison to both 
IShen et all i|200ot ) and IBaldrv et ail ||2012T ) is fully consis- 
tent, in that both studies also derive -R50 by fitting Sersic 
laws to the azimuthally- average surfa c e brig htness profiles 
of the observed galaxies. IReves et alj (|201ll ). on the other 
hand, perform somewhat more complex surface brightness 
modelling (they fit a 2D double exponential), which we do 
not attempt for the simulated galaxies. 

Over the range 9.0 < log 10 [M*(M Q )] < 10.4 the sim- 
ulated galaxies lie approx ima tely between the me dian rela- 
tions of lShen et all (|2003h and lBaldrv et all |2012f ) and have 
significant overlap with both (bottom pa nel) . The simulated 
galaxies are also consistent with data of I Reyes" et all (|201ll ) 
over this stellar mass range (taking into account error bars 
and scatter), although the trend with mass starts to deviate 
from the observations at masses of log lo [M*(M0)] ~ 10.0. 
The general agreement at low-moderate masses is encour- 
aging and, together with the agreement with the rotation 
curves shown above, indicates that the simulated galaxies 
do not suffer from the well-known 'angular momentum prob- 
lem'. This is not due to resolution, as results do not change 
significantly if we adopt the intermediate-resolution GIMIC 
simulations instead (see the Appendix). 

For log 10 [M*(MQ)] > 10.5 the simulated galaxies are 
obviously too compact compared to observed galaxies of the 
same stellar mass. This is consistent with our findings from 
the Tully-Fisher relation and the rotation curves in bins of 
stellar mass presented above. 

7 Naively, one might expect spheroid-dominated, low-angular 
momentum galaxies to have more peaked circular velocity curves 
than disc-dominated galaxies. This would indeed be true for 
galaxies of fixed stellar mass fraction. However, at fixed stellar 
mass the situation is reversed, since at the high mass end disc- 
dominated galaxies tend to have a higher stellar mass fraction 
than spheroid-dominated galaxie s (see Fig. 4), whic h is found to 
be the case observationally (e.g., Iputton et al. 2010). The higher 
stellar mass fraction leads to more peaked circular velocity curves 
for disc-dominated galaxies. 
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Figure 3. Comparison of simulated and observed relations be- 
tween galaxy stellar mass and size, as characterised by the radii 
which enclose 80% (top panel) and 50% (bottom panel) of the 
i-band flux in projection. The simulated galaxies lie approxi- 
m ately between the medians of the observed M 4 — R^q relations 
of Shc n et al. | (|2003h (fro m SDSS data , green curves with dashed 
lines indicating ±ler) and lBaldrv et al.l (2OI2I ) (from GAMA data, 
purple lines, with dashed lines indica ting dzlfr) and app roximately 
on top of the M* — Rgo relation of lReves et alj l|201lh over the 
range 9.0 < log 10 [M*(MQ)] < 10.4. At higher masses, the simu- 
lated galaxies are too compact compared to observed galaxies of 
the same stellar mass. 

3.4 Stellar mass - star formation efficiency 
relation 

We have shown that over the range of stellar mass, 9.0 < 
log 10 [M*(M Q )] < 10.5, the GIMIC high-resolution simula- 
tions approximately reproduce the observed TF relation as 
well as the rotation curves and galaxy sizes as a function 
of stellar mass. Given that for this mass range the rotation 
velocity in the simulations is dominated by dark matter, the 
agreement with the observed TF relation and rotation curves 
would imply that the models should have approximately 
the correct SF efficiency, esf = {M, /M200) / '(f2f,/f2 m ). Here 
M* is the total current stellar mass of the galaxy (which 
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Figure 4. Star formation efficiency, esFi as a function of galaxy 
stellar mass. Each filled circle represents a simulated galaxy, 
colour-coded by the fitted Sersic index, n. The dotted or- 
ange curve corre sponds to the (scatter- free) relation derived by 
iGuo et al.l ((2010) by abundance matching a stellar mass function 
derived from SDSS data with the halo mass function from the 
Millennium Simulation. The dashed orange curve is the result 
when a scatter of 0.15 dex in M* at fixed M200 is adopted in 
the abundance matching and egF is computed using the mean 
M200 in bins of stell ar mass. The hatche d black and red regions 
correspond to fits of iDutton et al.l l|2010h to 'direct' measures of 
the SF efficiency (from stacked weak lensing and satellite kine- 
matics measurements) of disc and elliptical galaxies, respectively 
(see text). The simulated disc galaxies match the observed egp of 
disc galaxies over the mass range 9.0 ^ log 10 [M*(MQ)] < 10.5. 
At higher masses, the simulated elliptical galaxies have too high 
£SF relative to what is inferred observationally. 



excludes surrounding satellites), M200 is the total mass 
(gas+stars+dark matter) within r2oo, and flb/^lm ~ 0.165 
is the universal ratio of baryonic to total matter densities 
(|Komatsu et al.ll201ll ). The SF efficiency is therefore the 
fraction of the total mass accounted for by stars, relative to 
the maximal case where all the baryons are converted into 
stars. Cosmological hydrodynamic simulations have had a 
notoriously difficult time in forming normal galaxies with 
the correct esf, generally generating values of e sf that are 
a fac t or of a few larger t han observed (see, e.g., Guo et all 
l20ld , IDutton et al1l2010l , and lAvila-Reese et alj 



20111 . who 



compare the results of a number of recent hydrodynamic 
simulations with observations). This is often referred to as 
the 'overcooling problem' of cosmological simulations. How 
do the GIMIC high-resolution simulations fair in this regard? 

In Fig. [4] we plot the trend between esf and stellar 
mass for the simulated galaxies at z = 0. Each of the 
dots in Fig. [4] represents a simulated galaxy. The dots are 
colour-coded according to the fitted Sersic index. Note that 
we include both disc- and spheroid-dominated galaxies in 
thi s plot. The dotte d orange curve is the relation derived 
bv lGuo et all (|201fj| ) by applying the abundance matching 



techniqu e to the SDSS g alaxy stellar mass function^ de- 
rived bv lLi fe White] ( 2009) and the dark matter halo mass 
function from the Millennium Simulation. This result as- 
sumes no scatter in the relation between M* and M200. The 
dashed orange curve is the result when a scatter of 0.15 
dex in M* at fixed M200 is adopted and esf is computed 
using the mean M200 in bins of stellar mass (Guo, priv. 
comm.). The sha ded black and red r egions correspond to 
the (by eye) fit of IDutton et al.1 (|2010l ) to 'direct' measure- 
ment^ of the SF efficiencies of disc and elliptical galaxies, 
respectively. The direct measurements are bas e d on stacked 
weak lensing l|Mandelbaum et al.|[200rj 120081; ISchulz et all 
20101) and stacked satellite kinematics l|Conrov et al.L 2007: 



Klypin fc Pradall2009l ; lMore et alfcoill) (hereafter WL/SK) 



analyses of disc and elliptical galaxies. By stacking the 
WL/SK of galaxies of fixed stellar mass (or luminosity), 
these authors derive the mean halo mass (e.g., M200) as a 
function of stellar mass and thereby co nstrain the SF effi- 
ciency. Importantly, IDutton et ail l|2010l ) have taken care to 
adjust the stellar masses from these studies so that they all 
correspond to the Chabrier IMF, and to e nsure that the halo 
mass is defined as M200 throughout (see IDutton et aill2010l 
for details). Note also that the stellar masses derived in the 
individual WL/S K studies are ty pic ally based on the singl e 
colour scalings of lBell et al.l i|2003l ) or lBell fc de Jond l|200ll ). 
We therefore apply a small mass-dependent adjustment so 
that stellar masses from these studies are consistent with 
those derived from five-band SDSS data l|Blanton fc Roweisl 
2007) (as described in Section 1). It is therefore possible 
to directly compare the WL/SK estimates of the SF effi- 
ciency (the shaded regions) with the abundance matching 
results represented by the dashed an d dott ed black curves. 
The shaded regions of IDutton et alj {2010) bracket the 2- 
sigma (95% CL) errors quoted in the individual WL/SK 
studies. 

Consistent with the results derived in Sections 3.1 and 
3.2, the simulated galaxies have SF efficiencies that are 
quite comparable to what is derived from stacked WL/SK 
analyses of disc galaxies over the range in galaxy mass, 
9.0 < log 10 [A'/*(M©)] < 10.5. The correct fraction is 
achieved through efficient (but energetically feasible) SN 
feed back which ejects gas from the system (see, e.g., Fig. 
7 of ICrain et ail |2010| ). This is ultimately what also al- 
lows the simulations to reproduce the TF relation and the 
shape and normalization of the rotation curves over this 
mass range. We point out that unlike some re c ent studies of 
indiv i dual galaxies (e.g ., IStinson et alj |2010| ; iMaccio et al.l 
|2012| ; iBrook et al] [2012). we have not had to invoke multi- 
ple sources of feedback in order to reproduce the observed 
SF efficiency. All that has been done is to adjust the mass- 
loading of the SN feedback so that the simulations roughly 



The stellar masses of Li fe White! d2009l) were derived using the 
five-band SDSS photome tric data (and redshift) as described in 
iBlanton fc Rowei 1 ll2007f) and assuming a Chabrier IMF, consis- 
tent with what is adopted in our cosmological simulations. 
9 We note the 'direct' measurements of halo mass are model- 
dependent, in that they assume a NFW halo profile with a par- 
ticular mass-concentration relation that does not account for pos- 
sible correlations between the scatter in concentration and galaxy 
morphology / colour . 
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reproduce the peak of the cosmic SF rate density of the 
Universe. 

Interestingly, most galaxies in the simulations fall below 
the relation derived from applying the abundance match- 
ing technique to the galaxy stellar mass fun ction (with or 
witho ut scatter included) for M* < 10 10 M^. iDutton et al.l 
l|2010h previously pointed out the factor of « 2 inconsis- 
tency between the abundance matching and direct measure 
studies. It is presently unclear what the origin of the discrep- 
ancy is. The discrepancy cannot be explained by differences 
in the cosmological para meters adopted f or the Millennium 
Simulation (from which iGuo et al] |2010| der ive their halo 
mass function) and current best estimates, as lDutton et al.l 
(2010) have shown that the discrepan cy exists for other 
abundance matching -based results (e.e.. iMoster et al"1l20ld : 
iBehroozi et aljfeoicf ) which adopted parameter values more 
in line with those derived from recent CMB analyses. Fur- 
thermore, the discrepancy is not simply due to the fact that 
the abundance matching technique is applied to the stellar 
mass function without regard to galaxy morphology, since 
disc galaxies dominate the stellar mass function at stellar 
masses where the discrepancy exists (thus, the abundance 
matching results should reflect those of disc galaxies at these 
masses). 

However, the fact that the simulations approximately 
match the TF relation and velocity rotation curves over 
the same range of stellar mass, suggests that the system- 
atic offset may lie mostly with the abundance matching 
results. Consistent with this hypothesis is the normalisa- 
tion offset between the TF relation deduced from the abun- 
dance matching t echniquq 10 ! and th e observed relation (see 
IGuo et al.ll2010h . IGuo et all (^OloT ) suggested that the off- 
set could be due to the fact that their calculations did not 
take into account adiabatic contraction of the dark matter 
halo due to centrally-concentrated cold baryons, but an al- 
ternative explanation is tha t the dark matter has not signif- 
icantly contracted (see, e .g. , iDuffv et al.ll2010l ; IDutton et al.l 
l201ll ; (Macci6 et al.ll2012i and Section 3.5 below) and that the 
abundance matching results systematically overestimate the 
stellar masses by a factor of ~ 2. 

Indeed, the agreement with the TF relation in the top 
panel of Fig. Q] implies that the simulated galaxies have the 
correct stellar mass fractions within r$o over the range 9.0 < 
log 10 [M„ (M Q )] < 10.5, since Vso is a direct measure of the 
total mass (stars+gas+DM) within that radius (under the 
assumption that non-circular motions are small). Assuming 
that the observed and simulated galaxies have similar (dark) 
matter distributions beyond this radius (which is implied 
by the WL/SK data since these shows that an NFW profile 
(N avarro et al.lll996l , ll997f i works well at these radii), there 
is therefore consistency between the match to the WL/SK 
data in Fig. [4] and the TF relation in the top panel of Fig. [I] 

A possible systematic error in the abundance matching 
formalism is the implicit assumption that satellite and cen- 



tral gal axies of fixed halo mass have the sam^| stellar mass. 
Indeed. [Rodrfguez-Puebla et alj l|201lf ) have shown that by 
excluding satellites from the abundance matching calcula- 
tion (i.e., both from the halo mass function and the observed 
stellar mass function) the star formation efficiency of central 
galaxies decreases slightly, which goes in the right direction 
to resolve the difference between the abundance matching 
results and those based on WL/SK. However, the effect ap- 
pears to be too small (only w 10 - 20%) to reconcile the 
difference. 

There is a small number of apparent outliers from the 
main esf — M* trend at low to moderate stellar masses. We 
have investigated the nature of these outliers and found that 
they always correspond to galaxies in close proximity to a 
massive galaxy group. The dark matter associated with the 
galaxy is clearly truncated (which is why esF is higher than 
typical). This may be because these galaxies have passed 
through the more massive neighbour (and hence have been 
tidally stripped) in the past, but have since come back out 
to larger radii (see the population of "backsplash" subhaloes 
on extreme orbits identified in cosm ological simulations by 
iGill et al.ll2005l ; lLudlow et al.ll2009l ). 

Returning to the general trend, for stellar masses 
log 10 [M„ (M©)] > 10.5 the simulations clearly have too 
high stellar mass fractions with respect to both the abun- 
dance matching results and the stacked measurements of 
disc and especially elliptical galaxies. Furthermore, the sim- 
ulations do not produce a bi-modal distribution of stel- 
lar mass fractions (corresponding to discs and ellipticals) 
over the range 10.5 < log lo [M„(M0)] ^ 11.0 as implied 
by stacked WL/SK measurements. It is conceivable that 
feedback from AGN, which was not included in our sim- 
ulations and which is expected to become relatively more 
important at high masses (and for galaxies with relatively 
large spheroidal components ) , may resolve these issues. In- 
deed, [McCaxthxllal] (|2010r i have used simulations from the 
Over whelmingly Large Simulations project |Schave et al.l 
2010), which uses the same simulation code as gimic, to 
show that the inclusion of AGN feedback in cosmological 
simulations can reduce the stellar mass of the most massive 
galaxies by up to an order of magnitude (see Fig. 7 of that 
study), yielding agreement with observations on the mass 
scales corresponding to groups. Finally, we point out that 
current state-of-the-art very high resolution zoomed simula- 
tions (that neglect AGN feedback) also suffer from overcool- 
ing at t hese high masseq 12 ! . as d emonstrated by the Aquila 
project (|Scannapieco et al]|201l| y 

3.5 Stellar mass - V80/V200 relation 

Related to the star formation efficiency, esf, is the ra- 
tio of galaxy rotation velocity to its virial circular ve- 
locity, which we define her e (following I Reyes et alj |2012| ; 
see also IDutton et all l2010t ) as Vgo/V2ooc, where V200C = 



IGuo et al. (2010) use abundance matching to derive a rela- 
tion between M, and M200 an d use the relation between V ma x 
and M200 from the Millennium Simulation to derive the relation 
between M* and V max - If dark matter dominates the rotation 
curves, then this is equivalent to the TF relation. 



11 Note, however, that usually abundance matching is performed 
for satellite galaxies at the time of accretion. 

12 The code used for the GIMIC simulations was also used in the 
Aquila project, which simulated a sing le M200 = 1-6 X 1O 12 M 
halo. The results of that individual galaxy are completely in line 
with the trends present in Figs. 1 and 3 of the present study. 
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Figure 5. Comparison of simulated and observed relations be- 
tween galaxy stellar mass and the ratio of the 'optical rota- 
tion velocity' (denned here as V$q) to virial circular velocity, 

V20OC = \J GM200 /)"200 > f° r disc galaxies at 2 0. The shaded 
regions encompass the ±1 sigm a (fine shading) an d ±2 sigma 
(spare shading) distributions of iReves et al. I j2012h . The simu- 
lated galaxies lie approximately on top of the observed relation 
over the wide range 9.0 < log lo [M*(M )] < 10.5. At higher 
masses, the simulated galaxies have circular velocities that are 
too large with respect to their virial circular velocities. 



y GM200 / 7*200 • If baryonic processes do not significantly 
affect the total matter distribution at small radii (i.e., so 
that it approximately traces an unmodified NFW profile) 
we would expec10 Vgo/p2ooc ~ l-max/Vaooc ~ 1.16 typically 
for this range of masses. If, on the other hand, baryons signif- 
icantly contract (expand) the potential well on scales similar 
to the optical extent of the galaxy, this will have the effe ct 
of raising (lowering) the K.pt/VW (e.g., |Puffv et al.ll2010l ). 

In Fig. [5] we examine the relation between galaxy stellar 
mass and V op t/Kdr = Vgo/Viooc for our simulate d disc galax- 
ies. W e make comparisons to the observations of lReves et al.l 
(2012), who measured this relation for a disc galaxy sam- 
ple that is repres entative of the Tully-Fishe r samp le of 
IReves et all l|201ll ) (see Fig. [TJ. IReves et all lj2012f ) em- 
ployed stacked galaxy-galaxy lensing to measure a mean 
M200 (and thus a mean V^ooc) in bins of stellar mass. 
Consistent with our previous results, the simulated galax- 
ies lie approximately on top of the observed relation for 
9.0 < log 10 [M»(M©)] < 10.5, while at higher masses the 
simulated galaxies have too high galaxy circular velocities 



13 Here we used the 2 = Millennium Simulation halo catalog 
to measure the median value of Vmax/V^oOc (1-16) for haloes in 
the mass range 11.0 < log 10 [M 2 oo(M Q )] < 13.0. Note that the 
ratio Vmax/V200c varies with halo mass, owing to the fact that 
lower mass haloes are more concentrated on average, but the de- 
pendence is very weak. 



with respect to the virial circular velocity of the halo in 
which they live. 

The red dashed line in Fig. [S] represents the ratio ex- 
pected in the absence of significant contraction or expan- 
sion of the inner potential well. If we restrict the com- 
parison to the mass range 9.0 < log 10 [M*(MQ)] < 10.5, 
both the GIMIC simulations and the observations are consis- 
tent with a mild degree of expansion at low stellar masses 
(log 10 [M*(Mo)] < 9.5) and a similarly mild amount of con- 
traction at higher stellar masses (log 10 [M*(Mo)] > 10.0). 
At very high masses (log 10 [M*(Mo)] > 10.5), the simulated 
potential wells have clearly undergone too much contraction 
as a result of overcooling. 

3.6 Stellar mass - sSFR relation 

As discussed in the Introduction, there has been consider- 
able progress in recent years in the ability of simulations 
to reproduce the structural and dynamical properties of 
normal disc galaxies at z = 0. This is due in large part 
to the implementation of strong feedback mechanisms in 
the simulations. However, it has been pointed out that in 
spite of this progress, there remain a number of impor- 
tant problems, including the apparent inability of simula- 
tions with strong feedback to reproduce the sSFRs of galax- 
ies locally an d at high redshi f t, particularly for low-m ass 
gal axies (e.g., Dave et al.ll201ll ; lAvila-Reese et al1l201ll ; but 



see iBrook et al.ll2012l ). Here we examine the trend between 
galaxy stellar mass and sSFR at z — for disc galaxies in 
the GIMIC simulations. 

In Fig. [S]we compare the relation between galaxy stellar 
mass and sSFR for disc g alaxies in the GIMI C simulations 
with the observations of ISalim et alj il2007l ). ISalim et al.l 
i|2007l) estimated star formation rates using UV data from 
GALEX for a large sample of ~ 50,000 optically-selected 
galaxies from the SDSS. The black shaded region in Fig. [5] 
represents the galaxy mass range for which the sSFRs are 
sensitive to numerical resolution and are likely not con- 
verged. 

At stellar masses of log 10 [M*(MQ)] > 10.0, the simu- 
lated galaxies roughly reproduce the observed trend, which 
is encouraging. Note that this (rough) agreement was not a 
given, as the feedback parameters in GIMIC have been tuned 
to match only the peak of the cosmic star formation rate 
density (at z ~ 1 — 3), not the star formation rates of galax- 
ies at the present day. Note, however, that the agreement at 
the highest masses (log 10 [M*(MQ)] > 10.5) must be fortu- 
itous, as we have already demonstrated that these galaxies 
suffer from overcooling (therefore the star formation rates 
of these galaxies are too high given their halo masses) . 

At lower stellar masses (log 10 [M»(M Q )] < 10.0) the 
simulated galaxie s fall below the observed relation of 
ISalim et alT(|2007l ). However, as we demonstrate in the Ap- 
pendix, the sSFRs of galaxies with such low masses are 
rather poorly converged at z = in our simulations. It is 
possible that a further increase in resolution would result 
in yet higher sSFRs, possibly resolving the discrepancy at 
low stell ar masses. Consiste nt with this hypothesis is the 
fact that lBrook et al' (2012) have used very high resolution 
zoomed simulations of low-mass galaxies and are able to re- 
produce the relation of lSalim et al.l (|2007l ) down to low stel- 
lar masses. Thus, it remains to be seen whether the problem 
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Figure 6. Comparison of simulated and observed relations be- 
tween stellar mass and the sSFR of disc galaxies at z = 0. The 
blue solid and dashed lines represent the median and ±1 sigma 
intrin sic scatter, respect ively, of the observed 'blue galaxy' sam- 
ple of lSalim et al.l (120071 1 . The black shaded region represents the 
galaxy mass range for which sSFRs are sensitive to resolution (i.e., 
are likely not converged - see the Appendix). For the simulated 
disc galaxies that are not currently forming stars (i.e., SFR= 0), 
we arbitrarily set their log 10 [sSFR]= —12.9 in this plot. There is 
rough agreement with the observed sSFRs for simulated galaxies 
with log 10 [M*(M o )] > 10.0. 



of some simulations in not reproducing the sSFRs of low- 
mass galaxies is a fundamental one. 



4 CONNECTING STAR FORMATION 
EFFICIENCY, THE STELLAR MASS 
FUNCTION, AND THE TF RELATION 

In this section we briefly comment on the link between the 
star formation efficiency (i.e., the M* — M*/M 2 qo relation), 
the galaxy stellar mass function, and the TF relation. 

The virial mass of a halo of mass M200 can be ex- 
pressed in terms of its virial circular velocity, V200 = 
(GM 200 /r 200 ) 1/2 , as 



M200 = 3.2 x lO^Mc 



V200 



100 kms~ 



(2) 



assu ming a Hubble constant of Ho = 72 km s 1 Mpc 
(e.g-. lMo et ai1ll998l ). Thus, if we knew the relation between 
galaxy stellar mass, M* , and M200 we could convert eqn. (2) 
into a TF relation. 

For halos with masses below M200 ~ 1O 14 M0, the halo 
mass function is well represented by a power-law of the form 



dn 
dM 20 o 



oc Af 9: 



(3) 



(e.g.. I Jenkins et al.ll200ll ; lReed et al]|2007l ). If we assume the 



galaxy stellar mass function also follows a power-law of the 
form 



dn 

1m2 



oc M~ a , 



(4) 



then matching abundances of haloes and galaxie^f] results 
in an esF-M2oo relation of 



M 2 



oc M, 



[(1.9-a)/0.9] 



(5) 



Inserting the above into eqn. (2) yields a TF relation of the 
form 



V200 oc Mi 1/3)(Q - 1)/a9 oc Mi Q " 1)/2 ' 7 



(6) 



Over the range 9.0 < log lo [M*(M )] < 10.0 (which 
avoids resolution issues at lower masses and the knee in the 
stellar mass function at higher masses, which would violate 
the power- law assumption in eqn. 4), the galaxy stellar mass 
function in GIMIC is characterised relatively well by a power- 
law with a ~ 1.6. From eqns. (5) and (6) this yields esf oc 
M^r 1 ^ 3 and V200 oc M^ ' 22 over the same range in stellar 
mass. These expectations match what we actually measure 
over this range in GIMIC well. 

It is therefore possible, at least in principle, to constrain 
the faint-end slope of the galaxy stellar mass function from 
the slope of the TF relation (or vice versa). On this basis, 
the observed TF relations in Fig. 1 suggest that a must be 
rel atively steep, > 1.5, which is close the value of 1.45 found 
bv lBaldrv et al.l 1 20121 ) for blue galaxies. 

A rigorous statement to this effect would necessitate re- 
laxing a number of assumptions. Firstly, we have assumed 
power-law dependencies for the halo and stellar mass func- 
tions over this mass range. This holds to a high level of accu- 
racy for the halo mass function but should be checked for the 
stellar mass function. Secondly, we have assumed that there 
is no scatter in the M*-M 2 oo relation. We have explored the 
effects of adding log-normal scatter to an idealised powerlaw 
relation between M» and M200 and find that it does not bias 
the recovered slope of the resulting stellar mass function. 
However, if the amplitude of the scatter were a strong func- 
tion of M200 and/or if the scatter were not log- normal, then 
this can result in a bias in a. Finally, we have derived the 
TF relation in eqn. 6 in terms of the virial circular velocity, 
whereas what is measured is closer to Vmax- If dark mat- 
ter dominates the dynamics, then to a good approximation 
Vmax ~ 1.2V200- This breaks down when applied to a very 
wide range of halo masses owing to the mass-concentration 
relation of dark matter haloes, an effect, however, which is 
unimportant for the mass range under consideration here. 
Potentially more relevant here is that baryonic physics can 
either contract or expand the dark ma tter distribution (e.g., 
iDuffv et al1l20ld ; iDutton et alll201ll ) and can therefore in- 
troduce a halo mass dependence in the relationship between 



14 Note that in Section 3.4 we argued there may be a systematic 
problem with galaxy stellar masses inferred through the abun- 
dance matching technique. However, this appears to be a problem 
in the normalisation, rather than in the slope (which is very simi- 
lar to that derived from stacked WL/SK). Here we are interested 
in how the slopes of the e — M* relation, the TF relation, and the 
stellar mass function are related. 
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Vmax and V2oo- In the GIMIC simulations, we find these ef- 
fects to be small over the mass range considered here (i.e., 
9.0 < Iog 10 [M*(M Q )] < 10.0). 



5 SUMMARY AND DISCUSSION 

We have used the GIMIC simulations to study the mass dis- 
tribution and SF efficiency of a large, cosmologically repre- 
sentative population of simulated disc galaxies over the wide 
range of stellar mass 9.0 ^ log lo [M*(M0)] ^ 11.5 at z = 0. 
These simulations include efficient, but energetically feasi- 
ble supernova feedback, but have not been tuned in any way 
to produce "realistic" disc galaxies. The main results of our 
study may be summarised as follows: 

• Over the stellar mass range 9.0 < log lo [A/„(M0)] < 
10.5, the simulations approximately reproduce: the observed 
Tully-Fisher relation (this is true for both the velocity mea- 
sured at the radius that encloses 80% of the galaxy's i- 
band light, Vgo(M*), and the peak velocity, Vmax(M,)), the 
observed rotation velocity curves of disc galaxies, the ob- 
served size-mass relation of disc galaxies, the relation stel- 
lar mass (M*) and optical rotation to virial circular veloc- 
ity ratio ('Kjpt/V^ir'), and the trend between SF efficiency 
and galaxy stellar mass, as inferred from stacked weak lens- 
ing and stacked satellite kinematics studies. Over this mass 
range, dark matter dominates the rotation curves of galax- 
ies. 

• A comparison of both the observed and simulated val- 
ues of 'Kjpt/Vvir' with that predicted by dark matter only 
simulations indicates that the potential wells of galaxies 



with relatively low stellar masses (log 10 [M, (Mg 



.5) 



have undergone a relatively mild amount of expansion, while 
the potential wells of galaxies with higher stellar masses 
(log 10 [M* (Mq)] > 10.0) have undergone a similarly mild 
degree of contraction (see Sec. 3.6). 

• For stellar masses log lo [M*(M0)] > 10.6, the simulated 
galaxies rotate faster than observed galaxies of the same 
stellar mass, are too compact, and have too large a fraction 
of their baryons locked up in stars (i.e., the simulations suffer 
from overcooling) . 

• The simulations also reproduce the specific star forma- 
tion rates of ~ L* galaxies (with log lo [M*(M0)] ~ 10.6) 
but have too low levels of SF in lower mass galaxies, which 
is plausibly due to the finite resolution of the simulations. 

• There is an intriguing discrepancy between the SF ef- 
ficiency of disc galaxies derived from abundance matching 
and 'direct' measurements derived from stacked weak lens- 
ing/satellite kinematics for M, < 1Q 1 O M0, which was pre- 
viously pointed out by iDutton et alj ([2010) . The fact that 
our simulations simultaneously reproduce the TF relation, 
rotation curves in bins of stellar mass, and the SF efficiency 
derived from direct measurements over this mass range, sug- 
gests that the unknown systematic may lie mostly with the 
abundance matching results. One possibility is that the im- 
plicit assumption that satellites (at accretion) and central 
galaxies of fixed halo mas s have the same star formati on 
efficiency is incorrect (e.g.. iRodrfguez-Puebla et al.ll201ll . 

• Simple scaling arguments presented in Section 4 (see 
also lMo et al.|[201Ql ) show that it is possible to constrain the 
faint-end slope of the galaxy stellar mass function from the 



slope of the TF relation (or vice versa). Current measure- 
ments of the TF relation suggest that the faint slope should 
be relatively steep (a > 1. 5), which is consiste nt with some 
recent measurements (e.g.. iBaldrv et al1l2012f ). 

The success of the GIMIC simulations in reproducing 
the observed trends at low to moderate stellar masses is 
directly att ributable to efficient SN feedback. As shown pre- 
viously by ICrain et all (|2010h . the SN feedback in GIMIC 
drives vigorous galactic winds, particularly at high red- 
shift, which eject a significant fraction of the baryons from 
haloes with masses log lo [M2oo(M0)] < 12.5 and pollute 
the intergalactic medium with metals. Thus, comparisons 
with observations of galactic outflows and the intergalac- 
tic medium represent promising further ways to test the 
model. We point out that efficient feedback was achieved 
without the need for extremely high resolution and is a re- 
sult of the way in which the feedback was implemented in 
the simulations. Recent zoomed simulations in which the 
feedback was inje cted thermally a nd th e cooling turn ed off 



temporarily (e.g., Governato e t al. 2010; Stinson et al ]l2010l: 



iGuedes et all l201ll ; iBrook et all [2012: Mac cio et all 120121 ) . 



appear to require much higher resolution before the feed- 
back becomes similarly efficient. 

For log lo [M*(M )] > 10.6, the galaxies in the GIMIC 
simulations suffer from strong overcooling, as do all other 
current cosmological simulations of disc galaxy formation 
that invoke only SN feedback. We conjecture that another 
energetically important form of feedback is required at high 



permassive black holes { 


Benson et al.[2003; Di Matteo et al. 


2005[:lBower et al.ll2006l; 


Croton et al.ll2006l:lBooth & Schave 


120091). iMcCarthv et al.l( 


201o|) have shown, using simulations 



that for m part of the Ove rwhelmingly Large Simulations 
project (jSchave et al.ll201fjh . that AGN feedback is crucial 
in reproducing the hot gas and stellar component of sys- 
tems with halo masses only slightly higher than explored 
here (corresponding to galaxy groups). It is therefore cer- 
tainly plausible that AGN are relevant over the mass scale 
where GIMIC (and other current simulations) fails. 

While the results of the present study certainly give 
reason for optimism, any successful model must not only 
reproduce these (and other) trends at z = but also their 
evolution to high redshift, which is by no means automatic. 
In a future study we intend to explore a wider range of 
observables (for both discs and ellipticals) over a wider range 
of redshifts and with a wider range of models, some of which 
include feedback from AGN. 
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Figure 7. Comparison of the Af* — M200 relations for the GIMIC 
intermediate- and high-resolution simulations. The filled circles 
correspond to the median halo mass, M200, in bins of galaxy 
stellar mass, M*, while the error bars correspond to the 14th and 
86th percentiles. There is reasonably good agreement when there 
are a minimum of ~ 100 star particles. We conservatively impose 
a cut of 500 star particles on the high-resolution simulation in the 
analyses presented in this paper. 
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APPENDIX: NUMERICAL CONVERGENCE 

Here we investigate the sensitivity of our results to the nu- 
merical resolution. In Fig. [7] we compare the M» — M200 
relations of the GIMIC intermediate- and high-resolution sim- 
ulations. In general there is good agreement when there is 
a minimum of ~ 100 star particles (modulo some sensitiv- 
ity to the exact mass scale where feedback stalls). We con- 
servatively impose a cut of 500 star particles on the high- 
resolution simulation (corresponding to log lo [M*(M0)] = 
9.0) in the paper. 

In Fig. [8] we compare the M t — Rso relations of the GIMIC 
intermediate- and high-resolution simulations. Again, there 
is reasonably good agreement when there is a minimum of 
~ 100 star particles, but with some differences occurring 
at the transition from the effective to ineffective feedback 
regimes. 

In Fig. [9] we compare the M„— sSFR relations of 
the GIMIC intermediate- and high-resolution simulations. 
There is reasonably good agreement for galaxies with 
log 10 [M»(M©)] > 9.7. At lower masses, the sSFR rate de- 
creases significantly with decreasing resolution. 
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Figure 8. Comparison of the M* — Rgo relations of disc- 
dominated galaxies for the GIMIC intermediate- and high- 
resolution simulations. The filled circles correspond to the me- 
dian value of RgQ in bins of galaxy stellar mass, M*, while the 
error bars correspond to the 14th and 86th percentiles. There is 
reasonably good agreement over the mass range considered here. 
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Figure 9. Comparison of the M*— sSFR relations of disc- 
dominated galaxies for the GIMIC intermediate- and high- 
resolution simulations. The filled circles correspond to the median 
sSFR in bins of galaxy stellar mass, M* , while the error bars cor- 
respond to the 14th and 86th percentiles. There is reasonably 
good agreement for galaxies with log 1( )[M*(Mo)] > 9.7. At lower 
masses, the sSFR rate decreases significantly with decreasing res- 
olution. 
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